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This work provides a thorough study of Levy, or heavy-tailed, random matrices (LMs). By 
analyzing the self-consistent equation on the probability distribution of the diagonal elements of 
the resolvent we establish the equation determining the localization transition and obtain the phase 
diagram. Using arguments based on supersymmetric field theory and Dyson Brownian motion we 
show that the eigenvalue statistics is the same one as of the Gaussian orthogonal ensemble in the 
whole delocalized phase and is Poisson in the localized phase. Our numerics confirm these findings, 
valid in the limit of infinitely large LMs, but also reveal that the characteristic scale governing finite 
size effects diverges much faster than a power law approaching the transition and is already very 
large far from it. This leads to a very wide crossover region in which the system looks as if it were 
in a mixed phase. Our results, together with the ones obtained previously, now provide a complete 
theory of Levy matrices. 


Since the well-known pioneering applications of Gaus¬ 
sian random matrices to nuclear spectra, random matrix 
theory (RMT) has found successful applications in many 
areas of physics [1] and also in other research fields such 
as wireless communications [2], financial risk [3], and bi¬ 
ology m- The reason for such remarkable versatility is 
that RMT provides universal results which are indepen¬ 
dent of the specific probability distribution of the random 
entries: only a few features that determine the univer¬ 
sality class matter. The most commonly studied RMs 
belong to the Gaussian ensembles pQ. They have been 
analyzed in great depth taking advantage of the symme¬ 
try under the orthogonal (or unitary/symplectic) group 
of the probability distribution. As an example of univer¬ 
sality, N x N real symmetric RMs, although they belong 
to the Gaussian orthogonal ensemble (GOE) only if the 
elements are Gaussian variables, display a GOE-like level 
statistics also when the distribution of the elements is not 
Gaussian, provided that it decreases fast enough to in¬ 
finity ehse 

There exists, however, a large set of matrices that fall 
out of the universality classes based on the Gaussian 
paradigm |7j. These are obtained when the entries are 
heavy-tailed i.i.d. random variables (i.e., with infinite 
variance). The reference case for this different universal¬ 
ity class corresponds to entries that are Levy distributed. 
This is the natural generalization of the Gaussian case 
since the limiting distribution of the sum of a large num¬ 
ber of heavy-tailed i.i.d. random variables is indeed a 
Levy distribution, as is the Gaussian distribution for 
nonheavy tailed random variables. Understanding the 
statistical spectral properties of these, so called, Levy 
matrices (LMs) is an exciting problem from the math¬ 
ematical and the physical sides [ME They represent 
a new (and very broad) universality class, with differ¬ 
ent and somehow unexpected properties with respect to 
the Gaussian case. Actually, a huge variety of distribu¬ 
tions in physics and in other disciplines exhibit power-law 


behavior. Accordingly, LMs appear in several contexts: 
in models of spin glasses with dipolar RKKY interac¬ 
tions EE in disordered electronic systems HE in port¬ 
folio optimization EE and in the study of correlations 
in big data sets [15] . just to cite a few. 

Contrary to the Gaussian case, the theory of random LMs 
is not yet well established. LMs were introduced in the 
pioneering work of Ref. Tj and further studied in Refs. E 
fl4l . By now, the behavior of their density of states is well 
understood (even rigorously) [7l-fT0j . Instead, on finer ob¬ 
servables, such as level and eigenfunction statistics, there 
are scarcer and even conflicting results. This is proba¬ 
bly due to the fact that the behavior of LMs is richer, 
and hence more difficult to understand, than the one of 
GOE matrices. For instance, a mobility edge separating 
high energy localized states from low energy extended 
states appears within their spectrum [TJ. It was also ar¬ 
gued that they display a new intermediate mixed phase, 
characterized by a nonuniversal level statistics. Although 
some aspects of the scenario put forward in Ref. [7] are 
in contradiction with recent rigorous results EE such 
mixed phase could indeed exist and actually be related 
to the one recently observed in the Anderson model on 
the Bethe lattice E3EE It may be the simplest case of 
the nonergodic delocalized phase advocated for quantum 
many body disordered systems in Ref. [21] . 

In the following we focus on N x N real symmetric 
matrices PL with entries hij = hji distributed indepen¬ 
dently according to a law, P(hij ) = N 1 ^ 
characterized by heavy tails: 

P ^ ~ 2N\Cj\ 1+fl ’ ^ °° 5 M < 2 ' 

The specific form of /( x) does not matter. For concrete¬ 
ness in numerical applications we will focus on a Student 
distribution with exponent 1 + /j, and symmetric entries, 
f(x) = f(—x). The scaling of the entries with N is such 
that almost all eigenvalues are 0(1) for N —> oo. 
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The first issue we address is determining the localization- 
delocalization transition line E*(p) in the E-p plane. 
In order to do so, we focus on the statistics of the di¬ 
agonal elements of the resolvent matrix G = \(E — 
irj)I — T-L}- 1 , which allows one to compute in the rj —► 
0+ limit spectral properties of T~L such as the global 
density of states p(E) = (1/iV) J2n=i — A„) = 
lim JJ _ > .o+ (1/Ntt) S Gu and the average inverse par¬ 
ticipation ratio (IPR) <T 2 , n > = (J2iLi l(*l n )| 4 ) = 

lim 7J _ >0 +(l/iV) J^iLi v\Gu\ 2 - As shown in Refs. [TUTU], 
the probability distribution Q(G) of a given Gu is ob¬ 
tained in the large -N limit from the equation: 

N 

G~^E-ir ] -Y. h % G n’ « 

i=i 

where all correlations between the terms on the rhs can 
be neglected and = denotes the equality in distribution 
between random variables. This leads to a self-consistent 
equation on Q(G), whose analysis yields the results on 
the density of states obtained in Refs. [3 HU]: For p<2, 
p(E) is a /x-dependent symmetric distribution with sup¬ 
port on the whole real axis and fat tails with exponent 
1 + p (the semicircle law is recovered for p > 2 only). 
There are several complementary ways to obtain the lo¬ 
calization transition from the statistics of the Gas. We 
have followed the one more likely to receive a rigorous 
treatment, as it was shown for the Anderson transition on 
the Bethe lattice [22]. It consists in studying the stability 
of the localized phase, checking whether adding a small 
imaginary part to Gu is an unstable perturbation 55] . 
Such stability is governed by an eigenvalue equation for 
the same integral operator found in Ref. jj], whose anal¬ 
ysis can be considerably simplified, as shown in Ref. [32j , 
and boils down to the following closed equation for the 
mobility edge E*(p), which is one of the main results of 
this work: 

K (si - a? /2 ) m*)\ 2 - 2Vv M(E*) + 1 = 0, (2) 

where = pT( 1/2 — p/2) 2 /2, = sin(np/2) and 

1(E) = J 0 +o ° e lkE dk/n. The function 

L C ^' > ' l3i ' E \k) is the Fourier transform of the probability 
distribution of the real part of the self-energy, that pre¬ 
vious works have shown to be a Levy stable distribution 
with exponent 1 + p/2 and parameters C(E ) and /3(E) 
determined self-consistently mmm- This equation has 
a solution for p £ (0,1) only. For p —>• 1 we find that 
E*(p) diverges as (1 — /x ) _1 . In Fig. [l] we show the nu¬ 
merical solution of Eq. Q for several values of p (we only 
consider E > 0 since the spectral properties are symmet¬ 
ric around zero). This quantitative phase diagram is in 
agreement with the sketch of Ref. [Tj and the numerics 
of Ref. [U] (except for p > 1 where the results were likely 
inaccurate due to the very large values of E that had to 



FIG. 1. Phase diagram of LMs in the p-E plane. 


be explored). 

We now address more subtle issues related to the level 
and eigenfunction statistics. We present first, two analyt¬ 
ical arguments which show that the statistics is a GOE in 
the whole delocalized phase and Poisson in the localized 
phase for N —» oo. The former is based on the super- 
symmetric zero-dimensional field theory introduced for 
random GOE matrices [23] ■ Since we follow closely the 
techniques developed in Refs. [25] [26] we just discuss the 
main steps and refer to Ref. [32, and a longer paper [27] 
for more details. The starting point is the field theory 
z = JUi cltlqe 5 ^!, with the action 


S = ^ M C ( E5 lm - h lm)<S>m + ^ <f>l 

l,m l 


- i0+ 


N 


The field f I» ?: is a eight-component super-vector 

($f ) ,$f ) ) = (S?,Slxi,ti,P?,P?,rk,rii), where each 

( 12 ) 

of the four component supervector $) is formed by 
two real and two Grassman variables. The matrix £ is 
diagonal with elements (1,1,1,1, —1, —1, —1, —1). The 
level statistics, in particular the density of states and 
the correlation between two levels at distance r/2N, can 
be obtained from correlation functions of the fields [54] . 
Averaging over the the matrix elements and introducing 
the function p($) = — $*) one can rewrite Z 

as fVp(^>)e s ^ with the action reading 


S = 


-NE 


1 / 

+ 2 (r 


— JV fd$p($) log/)(<!*) + ) NJ< d t I*d'I'g($)C( f l>l£'I')p('I'), 


where C(y) = pj 2 |4i+m [ ex p(— ixy) — 1]. Since the sec¬ 
ond term is subleading compared to the other three that 
are O(N), one can neglect it at first and perform a saddle 
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point. The solution of the corresponding equation reads 

/j($) = J dEi?(E) exp 

where, as it can be shown in full generality [25} [32], i?(E) 
is the probability distribution of the local self-energy, 
which coincides with the complex Levy stable law rig¬ 
orously proven in Ref. m (see Ref. (32]). Note that the 
saddle point equation is invariant under the symmetry 
<f> —► T'L where the super matrix T verifies the equation 
T^CT = 1. Thus given a solution p($), Pr(^) = p(T3>) 
is also a solution. The localization transition corresponds 
to the breaking of this symmetry [23} [25]: in the local¬ 
ized phase the typical value of the imaginary part of the 
self-energy is zero, whereas it is finite in the delocalized 
phase. In consequence, in the former case p( < f > ) is a func¬ 
tion of $!£<!> only, invariant under the symmetry gener¬ 
ated by T, whereas in the latter it depends also on $1$. 
Since this dependence breaks the symmetry there is a 
manifold of solutions p-r(^). It is the integration over 
this manifold that leads to GOE statistics for the level 
correlations. The derivation is identical to the one pre¬ 
sented in Ref. [25] since the only term in the action that 
depends on 7”, i.e., that breaks the symmetry, is the r 
one as it happens for Erdos-Renyi graphs [28] and GOE 
RMs |24j . In the localized phase, the saddle point solu¬ 
tion is instead unique. Therefore no integration over T 
has to be performed and this leads to uncorrelated levels, 
i.e., Poisson statistics [251 . 

Let us now turn to the other analytical argument, which 
is very straightforward but limited to /i > 1 only. Taking 
inspiration from the recent mathematical breakthrough 
on RMT (5], we slightly modify the distribution P(hij) 
into (1 - e)P(hij) + where W(x) is a 

Gaussian distribution with unit variance. This is equiva¬ 
lent to modifying TL into TL e = (1 — e)TL + eW where TL is 
a LM and W a very small GOE matrix whose elements 
have exactly the same scaling with TV than the ones of TL. 
Since this change does not alter the fat tails of the matrix 
elements, one naturally expects TL e and TL to be in the 
same universality class for any e < 1 and in particular 
for e —* 0. The statistics of the modified LM — and, by 
the previous argument, of TL —can be obtained using the 
Dyson Brownian motion (DBM): TLe can be interpreted, 
in the basis that diagonalizes TL, as a diagonal matrix to 
which an infinite number of infinitesimal GOE matrices 
have been added. The probability of the eigenvalues of 
TL e is therefore given by the DBM starting from the eigen¬ 
values of TL , and evolving over a fictive time of the order 
N~ 1 ^. Recent rigorous results [5] guarantee that the 
DBM has enough “time” to reach its stationary distri¬ 
bution, which is the GOE distribution, if TV -1 /^ TV -1 
and the typical level spacing of TL is 0(l/N) — a very 
reasonable assumption that agrees well with the numer¬ 
ics. This implies that for p > 1 the level statistics of 
the modified LM, and hence of the original LM too, is 



FIG. 2. ln[((r) — (r)p)/((r)GOE — (r)p)] as a function of 

E for different system sizes and for p = 1.5 (top) and p = 
0.5 (bottom). The dashed line represents the position of the 
mobility edge, E* ~ 3.85. 


indeed GOE-like in the bulk of the spectrum [31J (see 
Refs. [27} 22] for more details). 

We now present several numerical results with the aim 
of backing up our previous analytical arguments and also 
of studying the behavior of large but finite LMs. In 
applications TV is never truly infinite, actually in sev¬ 
eral cases it can be just a few thousand. Thus, it is 
of paramount importance to study finite size effects and 
determining the characteristic value of TV above which 
the TV —> oo limit is recovered. We performed exact di- 
agonalization of LMs for several system sizes TV = 2 n , 
from n = 8 to n = 15 and averaging over 2 22_rl real¬ 
izations of the disorder. We have resolved the energy 
spectrum in 64 small intervals u, centered around the 
energies E v = (A ra ) ng „, and analyzed the statistics of 
eigenvalues and eigenfunctions in each one of them. We 
have focused on several observables that display differ¬ 
ent universal behaviors in the GOE and Poisson regimes: 
The first probe, introduced in Ref. [25], is the ratio of 
adjacent gaps r n = min{5 n , 6 n +i}/ max{5 n ,<5 ra -t-i} where 
8 n = A n+ i — X n > 0 denotes the level spacings be¬ 
tween neighboring eigenvalues. It has different univer¬ 
sal distributions in the GOE and Poisson cases encod¬ 
ing, respectively, the repulsion or the independence of 
levels. The second one is the overlap between eigen¬ 
vectors corresponding to subsequent eigenvalues, defined 
as q n = Yh=\ |(*|n)||(i|n+ 1)|. Its typical value = 
gO'i r;„)„<?,. a n ows us t 0 m ake the difference between the 
localized phase, in which subsequent eigenvectors do not 
overlap (g typ = 0), and the delocalized GOE one in 
which they do (< 7 typ = 2/7r). Finally, the wave func¬ 
tion support set, recently introduced in Ref. [35]) is de¬ 
fined for an eigenvector n with sites ordered according 
to |(z|n)| > |(i + l|n)| as the sets of sites i < such 
that ]Tfli |(T|n)| 2 < 1 — e < ]T)f=i +1 l(*l n )| 2 - Tlie sca l“ 


- 3?E) + j 
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FIG. 3. Main panel: log N^(E) = -log (3Gg p p(£)) = 
— (log QGu) — log((AG)/7r) as a function of E for p = 0.5. 
Inset: q typ as a function of N for different energies and p = 
0.5, showing the position of N rn (E). 


ing of Se(N) for TV —> oo and e arbitrary small but 
finite allows us to discriminate between a localized and 
extended phase. The analysis of all these probes clearly 
shows that for p > 1 the eigenvalues and eigenvectors 
statistics is a GOE in the limit of large TV, in agreement 
with our previous arguments (and also with the rigorous 
results on the delocalized nature of wave functions H3I). 
As an example we show in the top panel of Fig. [2] the 
behavior of the average value of r n for p = 1.5. It clearly 
converges in the limit TV —> oo and for any energy E to 
the value (r n )coE — 0.53 characteristic of GOE statis¬ 
tics. All other probes show a similar convergence to the 
values expected for GOE statistics (see Ref. [32] for the 
corresponding plots). This is no longer true for p < 1, 
where the situation is more involved. In the lower panel 
of Fig. [2] we show again the behavior of the average of 
r n but now for p = 0.5. For small and large energies 
we find the values (r )goe — 0.53 and ( r)p ~ 0.39 cor¬ 
responding respectively to GOE and Poisson statistics. 
Moreover, the curves corresponding to different values of 
TV seem to cross much before the localization transition, 
that our previous analytical results located at E* ~ 3.85 
for p = 0.5. If this were representative of the truly 
asymptotic large- TV behavior then it would possibly sig¬ 
nal the existence of a mixed phase which could be de¬ 
localized but nonergodic , i.e., not displaying GOE statis¬ 
tics. However, analyzing carefully the data—thanks to 
the large number of samples used to average over the 
disorder—we find that the crossing point is in fact very 
slowly drifting towards higher energies as TV is increased. 
The same behavior is found for all the probes we stud¬ 
ied. As an example, in the inset of Fig. [ 3 ] we plot g* yp 
as a function of the system size for energies belonging 
to the crossing region. This indeed shows that < 7 tyP is a 
nonmonotonic function of TV. We can then define a char¬ 


acteristic matrix size, N m {E ), such that for TV -C N m (E) 
the statistics appears to be intermediate between Pois¬ 
son and GOE (see Ref. [32]), whereas for TV N m (E) it 
tends again toward GOE. 

The existence of a crossover size can be understood from 
the properties of the distribution Q(G). What character¬ 
izes the delocalized phase is that, at any site i, the imag¬ 
inary part of Gu receives an infinitesimal contribution 
from an infinite number of eigenfunctions. This leads to a 
typical value of Gp (defined as 3G*f p = e < lo s ^ w hich 
is finite for TV — > 00 and rj —> 0. Instead, 3G*f p = 0 in 
the localized phase. Approaching the transition from the 
delocalized side, 3G-f p becomes extremely small. Thus, 
one needs to take large enough systems in order to realize 
that it is different from zero, and hence that the system 
is in the delocalized and GOE-like phase. The argument, 
which is based on the interpretation of 3Gp as the local 
density of states, is as follows. The number of states per 
unit of energy close to E is Np(E). This number, mul¬ 
tiplied by the typical value of the local density of states, 
has to be larger than one in order to be in a regime repre¬ 
sentative of the large-TV limit. This defines the crossover 
scale N' m (E) oc l/(SG[f p p(E)). We have compared nu¬ 
merically \nN' m {E) and In N m (E) and found that they 
are indeed proportional (see Ref. [32] for a plot), thus 
showing that our argument correctly captures the origin 
of the finite size effects. We plot the crossover scale [ac¬ 
tually N' m (E)] as a function of E in Fig. [ 3 ] for p = 0.5: it 
diverges very fast approaching E*(p). A good fit is pro¬ 
vided by an essential singularity. These results therefore 
unveil what is the mechanism responsible for the non- 
GOE statistics observed for finite LMs in a wide regime 
before the localization transition. 

In conclusion, we have presented a thorough analysis 
of the eigenvalues and eigenvectors statistics of random 
Levy matrices. We have shown that the localization and 
the level statistics transitions coincide but also unveil the 
existence of a crossover scale which is very large even far 
from the transition. Thus, many practical cases are ex¬ 
pected to be in the TV <C N m (E) regime. In consequence, 
the mixed behavior proposed in Ref. IT] will be often 
present in practice even though it is absent in the large- 
TV limit. Our work, together with the results obtained 
previously, now provides a complete theory of LMs. 
There are several directions worth pursuing more. It 
would be interesting to determine analytically the form 
of the divergence of N m (E). On the basis of our numer¬ 
ics and in analogy with previous works [24] [26} [28] we 
expect N m {E) oc e c K E *~ E l a . Most probably, the emer¬ 
gence of the crossover scale producing an apparent mixed 
phase takes place in several other related situations (e.g., 
Ref. [19]) that are, therefore, to be reanalyzed. Finally, 
our results provide a guideline for mathematicians work¬ 
ing on RMT. Thanks to the recent advances in the math¬ 
ematical analysis of random matrices [5] and localization 
phenomena [22] our findings are likely to be rigorously 
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proven in a not too distant future. 
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EPAPS 


In this supplementary material we provide more details and results related to several points discussed in the main 
text. For definiteness we shall consider N x N real symmetric (Wigner) Levy Matrices H with entries hjj = 
distributed independently according to a student distribution with exponent 1 + ji and typical value of order N~^^\ 


p{h ij ) = e(\h ij \>N- 1 ^) 




2 N\hi 


I l+M 


(3) 


COMPUTATION OF THE MOBILITY EDGE 


The linearized recursive equations 


In order to determine the position of the mobility edge, it is convenient to introduce the self-energy = E — ir\ — 
G~3 = Si + iAj, and linearize the recursive equations for the diagonal elements of the resolvent matrix [Eq. (1) of the 
main text] with respect to their imaginary part in the limit 77 —>• 0 + : 


N 


N 


hUE-Sj) 


N 


h 2 

^ \ ' 'hj 


•5. - E h l -T. i E-s,r+( V +Ar 


. ,E-Sj 

7 = 1 J 


h 2 . 

,l i3 


N 


-E 


2 A U 


(4a) 

(4b) 


The recursion relation (20 1 for the real part of the self-energy totally decouples from the equation (4b) for its imaginary 


part. The marginal probability distribution of the real part can thus be computed using the generalized central limit 


theorem. One possible route to obtain the recursion relations [Eq. (1) of the main text and Eqs. (20) and (4b)] is 
provided by the cavity method [I]. Within this framework, the Gjj s appearing in Eq. (20) are the diagonal element 


of the resolvent of a LM of (TV — 1) x ( N — 1) elements, where the i-th row and column have been removed. As a 
consequence in the N — > 00 limit the matrix elements hij and the terms Gjj are independent and uncorrelated by 
construction. Since the variance of the entries is infinite for ji < 2, by virtue of the generalized central limit theorem 
the probability distribution of the variable S tends, for N —» 00, to a Levy stable distribution, L G j E ^’^ E \S), with 
stability index /z/2, and “effective range” C(E) and “asymmetry parameter” /3(E) given by: 


c=r ( i -f) c ° s (f) 


7=1 


0 = 


^EE sign(sRG^) 


(5) 


fEE|^. 


77 


In the N —> 00 limit the diagonal elements of the resolvent become independent and identically distributed. We can 
therefore replace the sums in the r.h.s. of Eqs. © by integrals over the marginal probability distribution QnffiG). 
In the limit of vanishing imaginary part, one has that Si = E — l/JiG^. Hence by changing variable one has: 


Qfl(3?G) d^G = (e - E) 


d!RG 

|5RG | 2 


As a consequence, Eq. © can be written as a set of two coupled self-consistent equation for the parameters C(E) 
and /3(E) [TJ(2], which can be solved numerically with a reasonably high degree of accuracy. 


C(E) = r (1 - |) cos (^) j^L G ^ E \S) \E~S\~^ dS, 
EE L^ )ME) (S) sign (E - S) IE - S\S dS 


( 6 a) 


P(E) = 


+00 t C(E),0(E). 


f 

J — 00 fj,/ 2 


\E-S\~*dS 


( 6 b) 

















7 


The mapping to directed polymers in random media 


We now focus on the recursive equation (4b I for the imaginary part of the self-energy, 
appearing in the sum of the r.h.s. by their expression in terms of their “neighbors 
obtains: 


If one replaces the AjS 
iteratively, say R times, one 


A _ V” 'Sia y' 'b 2 i 3 y' "'iR-iiR . 

11 ,k^ E - s ^k^ E - s ^''\ R k-A E - s ^ 2 iR ' 

As a result, the equation for the imaginary part of the self-energy can be interpreted as a sum over directed paths of 
length R originating from the site i\. For each edge (i n +i, z n ) crossed by a given path, the contribution to A^ coming 
from the path picks the random factor h 2 ni /(E — S in+1 ) 2 . The equation above can be then telescoped as: 


Ah 


E 


n 

(<«+1,*«) 


(E-S in+1 )*J 


A i R , 


(7) 


The imaginary part of the self-energy thus satisfies the exact same recursive equation as the partition function 
of directed polymers in random media (DPRM) [31 in presence of quenched (correlated) bond disorder = 

hfj/(E — Sj ) 2 . One can therefore use the results established in this context to analyze the properties of the distribution 
of and study the localization transition of LMs, as detailed below. 

It is well known that, depending on the strength of the disorder, DPs have a freezing transition (called one-step 
Replica Symmetry Breaking) akin to the glass phase of the Random Energy Model [3J :J] . The sum in 0 is over 
an exponential number of paths, (N — 1 )\/(N — R — 1)! ~ N R . In the large R limit, two cases are possible: the 
sum is dominated either by few paths that give a 0 ( 1 ) contribution, or by an exponential number of paths, each 
of them giving a very small contribution but such that their sum is 0(1). The freezing-glass transition of DPRM 
corresponds to the transition between these two regimes. It is then clear that in fact, mutatis mutandis, the glass 
phase of DPRM corresponds to the Anderson localized regime, where the imaginary part of the self-energy is of order 
77 with probability 1 on (almost) all the sites, except on extremely rare and distant resonances where it is of 0 ( 1 ) and 
the sum in ([ 7 ]) is dominated by very few paths. Conversely, the ergodic phase of DPs corresponds to the delocalized 
regime, where Aj is finite on all the sites and all paths give a non-zero contribution to the sum. 

Such transition is related to an ergodicity breaking. In order to determine the transition point, let us imagine to 
introduce n replicas of the system with the same realization of the disorder (i.e., n identical LMs with the same matrix 
elements). Following the analogy with DPRM, the “partition function” of the replicated system is 


Z n 


n 


e n 

V (*„+!,»„} 


and (minus) the (quenched) “free energy” per site is <f> = log Z n /Rn. (In the following for simplicity we will instead 
compute the “annealed” free energy <f> = log Z n /Rn, since it yields the same transition point as the quenched one). In 
the glass phase the sum in Q is dominated by few, 0(1), paths. One can thus assume that the n replicas are divided 
in n/m groups of m replicas all freezed in the same specific path (one-step Replica Symmetry Breaking ansatz). One 
than has: 


Z n = 


e n e 


In this case the (annealed) replicated free energy of DPRM reads: 


, * (ra - E) = is; l0S (£ n 

V (*n+l»*n 




E-Sd 


i n +i 


( 8 ) 


which gives the behavior of the typical value of A” 1 . The free energy needs then to be extremized with respect to 
the parameter m: d(f>/dm\ m=rn * = 0. If one finds that cf>(rn*,E) < 0, then this implies that the partition function is 
















exponentially small and that the typical value of the imaginary part of the self-energy vanishes exponentially under 
iteration. Hence the system is localized. Conversely, for <f)(m*,E) > 0 the typical value of the imaginary part of the 
self-energy grows exponentially under iteration and the system is in the delocalized phase. The Anderson localization 
is therefore given by: 


T d(j){m , E ) 

< dm m— 

{ <t>(m*,E) = 0. 


For the Anderson model on the Bethe lattice the localization transition takes place for m = 1/2, as it has been 
rigorously proven in El El and indirectly found in [6]. It turns out that the this is also the case for Levy Matrices. 


The exact equation for the mobility edge 


Considering the site i n of a given path, using the recursion relation (20), one can rewrite the real part of the 
self-energy as: 




N-2 ;2 

E-S' 


i'=1 


where the sum over i' n runs over all the N — 2 neighbors of i n except the sites i n and i n +\ which belong to the path. 
As a result, averaging over the quenched disorder the free energy ([8]) is written in terms of the largest eigenvalue, 
A(m, E), of the following transfer-matrix integral operator: 


N-2 


Z in (SiJ = (N- 2) dS in+1 Z in+1 (S in+1 )dh in , in+l P(h in>in+1 ) [] [dh in>i>n P(h ln yJ dSi' n L°WM E \Si'J 


n n ,ln+l 


E ~ S in+1 


l h i i +1 ^ h - 


(9) 


E-S* 


where the factor N — 2 accounts for the number of way one can choose the neighbors i n +i among all the N — 1 
neighbors of i n except the site i n -i- For large R one has that <j>(m,E ) ~ (1/to) log A(m, E). As a consequence, the 
mobility edge is found at the value of E* where: 


1 

TO 


log A(to, E*) = 0, 


d_ 

dm 


1 

to 


log A (to., E*) 


= 0 . 


This yields: 


r A(to,£*) = 1, 

{ d (io) 

£*)=». 

Note that this is equivalent to study the stability of the localized phase by checking whether a small imaginary part 
vanishes exponentially under iteration, as done in (in this case A (to, E) corresponds to the Lyapunov exponent 

of the imaginary part of the self-energy and to plays the role of the exponent of the power-law tails of its marginal 
probability distribution). It is convenient to introduce the variable S via the relation 
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In the thermodynamic limit, N — 2 ~ IV — 1 ~ IV, according to the generalized central limit theorem S has the same 
stationary distribution of the real part of the self-energy, [S)\ 


. N-2 


N ~ 2 h 2 


i' n y *ii=i lr 

Eq. § then becomes: 

Z(X) = N j dhP(h) dSL^f )ME) (S) dX'Z(X') 5 (x-S 


N—too W 2 


h 2 


E-X’ 


E-X’ 


2m 


In order to solve the eigenvalue problem, we take the Fourier transform of both sides of Eq. (Ill and obtain: 

2m 

Z(X')e- ikh2/{E - x,) , 


Z(k) = NL E j E)ME \k) J dhdX' P{h) 


E-X' 


where L E j E ^’^^ E \k) is the Fourier transform of the Levy stable distribution: 


Using the fact that: 


C(E),p(E)^ = exp _c(E)\k\^ /2 (l + 0(E) tan sign(&)) 


poo ikx 

dx -= e i 5( 1 -“) sisn ( fe )|A;|a-i r ( 1 - a). 


( 11 ) 


( 12 ) 


(13) 


(14) 


we can easily integrate over the disorder matrix element distribution obtaining: 

r oo 

/ dhP{h) \h\ 2m e ~ ikh2 /( E ~ x ') = JJx r(m — m/ 2) e -i f (m-/i/2 ) sign(fe ( B-x,) ) 
Jo 2IV 


E-X' 


ti/2-r 


Plugging the result above into Eq. (12) we get: 


/*+ o° 

m = f r(m - At /2) \kr /2 ~ m L E j E) ^ E \k) / dX' 

" J — oo 


Z(X') 


\E-X'\ m +»/ 2 


g-»f (m-/i/2)sign(fe(B-X')) 


(15) 


Note that this is exactly the same equation found in [T| , where the authors studied the stability of the localized phase 
by determining the Lyapunov exponent of the imaginary part of the self-energy under iteration. 

It is convenient to replace Z(X') by the inverse Fourier transform of Z(k'). We then perform the integral over dX' 
by separating it into two pieces, and changing variable E — S —> z in the interval (—oo ,E) and S — E —> z in the 
interval (E,+ oo). Using again Eq. (Tbfl) we find: 


p+OO 


dX’ 


Ak'X' 


E - X'\ m +M/2 


e -if (m-/i/2)sign(fc(B-X')) _ 


e ik 'E _ m _ fj^/2) IJfe 7 1 T "-+i Lt / 2 — 1 g-*f (m->i/2)sign(fe) g -if (l-m-p/2)sign(fc') gif (m-/j/2)sign(fe) (l-m-/t/2)sign(fc') 


We plug back this last result into Eq. ( fT5[ ) and get: 

Z+{k) = f r(m - m/2) r(l - m - m/ 2) |fcr /2 " m ^ (E) (fc) 


#t/2 

Z_{k) = ^ T(m — m/2) P(1 — m — m/2) L E j E ' > ’ /3(E ' 1 (k) sin (7 tto) 1+ + sin I- 


sin -1+ + sin (7rm) J_ 


/7TMX 


r+°° dh' , 

/+= / — e ifc E lUr+^Z 2 - 1 Z(Jk!), 

Jo 7r 

/ o J I,/ 

- gife E | fc ,| m+ ^/2-l _ 

7T 


where /+ and /_ are defined as: 
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We also have defined Z + {k) and Z_(k) as the function Z(k) restricted to the regions k > 0 and k < 0 respectively, 
and introduce the coefficients: 


K m ,, = ^ T ( m - 




) r ( 


1 — to — 




„ _ . / 7r A i \ 

S/i m V 2 / ’ 
s m = sin ( 7 rm) , 

as done in the main text. We multiply both Z + (k) and Z_(k) by e lkE |/ c | m +/ i / 2 - 1 and integrate them over dfc/-7r in 
the intervals (0,+ 00 ) and (— 00 ,0) respectively. We then obtain: 


1 + — A/n./i ( S 11 1 -\r T 


^m,/i 


+ ) 


I— — E m 11 (Sm-E|- T S^I—') ^— 


(16) 


where £+ and are defined as: 


/■-{-OO 


d k 


l + — 


i ._ = 


ikE 1 f C(E),f)(E) 




(k)- 


dk JkE,un-l fC(E),P(E) 


— e lkE I £f- 
7 r 


J fi/2 


(k). 


The 2x2 linear system (16) only has non-trivial solutions different from zero if the determinant of the matrix of the 
coefficients vanishes. Hence the equation (101 for the mobility edge reads: 


K 2 m J + t- [si -s 2 m ]- K m (£+ + iJ) S[l + 1 = 0. 


(17) 


a . s ~</ zp\ Of 

Using the specific form of L L '(&), Eq. (13), it is straightforward to show that i- = i* + . Interestingly enough, 
we remark that, as for the Anderson model on the Bethe lattice [6j, the l.h.s. of Eq. 0 is a symmetric function of m 
around m = 1/2. This implies that, if a solution of Eq. ( fl7| ) exists, then the stationary condition dX(m, E*)/dm = 0 
can only be verihed for m = 1/2. Hence, Eq. 0, finally becomes Eq. (2) of the main text. 

We have solved Eq. 0 — together with Eqs. ( [6a] ) and Eqs. ( |6b[ ) -numerically for several values of fL £ (0,1) and 
E (and for to = 1/2), and obtained the phase diagram of fig. 1 of the main text. 


SUPER-SYMMETRIC FORMALISM 

We give here more details on the super-symmetric formalism we discussed in the main text. The super-symmetric 
method to study random matrices is well established by now [S]. Moreover, our derivations are along the lines of the 
ones developed by Mirlin and Fyodorov for the connectivity matrix and the Anderson localization of finite connectivity 
random graphs [9/ For completeness, we present the main ideas and technical steps. A full derivation will be presented 
elsewhere m- 


Action in terms of p(<!>) 


As stated in the main text, the starting point is the field theory Z = /a dd^H with the action gj: 


5 = l - ( Y, - hlm)$ m + 

\ l,m l 

The field is a eight-component super-vector ($- 1 \$- 2 ^) = (£“, S'/, Xi, X*, -P“, Pi, ??*)> where each of the four 

( 12 ). 

component super-vector dq is formed by two real and two Grassman variables. The matrix C is diagonal with 
elements (1,1,1,1, —1, —1, —1, —1). In order to average over the distribution of the matrix elements one has to 
compute: 



e -ih lm <S>\C<S>. 

l<m 


I 
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Note that we have used explicitly that Using that 


-ih lm , A L _ f dhim / 

- + ivi 2 |^ to ii+m l e ; ■ 


2|^m| 1+ ^ 

one finds an action which reads (up to subleasing terms in TV, that do not play any role for our derivations): 
c 1 I r i ?~ + *u ' \ 1* \ -> / d h 


r + *o+\ , ft v- 
‘' Ar ' + 2TV ^. 

lm 


TV 


„-ih§t 


2|h| 1 +M 


- 1 


)• 


Introducing the function p(<E>) = [TV <5(d> — dq) one can rewrite 5 a as 


5 a = l -NE j d$p($)&£$ + l -(r + T0+) 


+ -NJ d$d'kp($)C , ($ t £'I')p(’I'). 


the function C(y) being the one defined in the main text: C(y) = p J 2 |^jf +M [exp {—ixy) — 1]. It still remains to 
integrate over all <3>iS. Since the action depends on the dqs through p(d>) only, one can first integrate over all d>,s that 
correspond to the same p(d>). This leads to an additional entropic-like term in the action —TV J dd>p(d>) logp($) (the 
computation is standard even though generically is not done with super-field). The final result is that the field theory 
has been now transformed in a new one: Z = J Dp(<£>)e s ^ where the integral is over all normalized p(<I>) with the 
action reading 


S[ P } = 


1 

2 


NE 



-TV 


dd>p(d>) logp(d>) + 



• (18) 


The interest of this formulation is that because of the TV which can be factored out in the action, one can evaluate it 
by the saddle-point method. 


Relationship between p(d?) and R{ E) 


Before discussing the corresponding equation it is useful to remark that the average value of /o(<E>), which corresponds 
to the saddle point of the previous integral, has a particularly illuminating expression in terms of the distribution of 
the local self-energy i?(E) [TT] . Before that the average over the disorder is performed, the field theory is Gaussian. 
Hence, by integrating all fields but T, one remains with a Gaussian integral to handle. Using that the field theory is 
constructed in such a way that ($- 1) l<I>( 1 ^) = 4 iGij and (d>- 2 ^d>^) = 4 iG*j (Gij is the resolvent and (•) denotes the 
average over the held theory at fixed disorder) the average (<S($ — dq)) turns out to be the Gaussian measure on d^. 
By collecting all terms and averaging over the disorder one finds: 

jfW) = ^ E ex p (^ f £<*>(£ - he«) + . 

i x ' 

By introducing the distribution of the local self-energy R( E), one gets the expression quoted in the text and already 
derived in HU: 

(p(T)) = J dEii(E) exp £$(£?-RE)+ . (19) 

Since the held theory can be solved by the saddle point method, the saddle-point value of p(d>) has to satisfy the 
previous equation as stated in the text. We show below that this is indeed the case. 


Equation on T?(E) 


As discussed in the main text, as well as in the hrst section of the SM, it was shown in [U H ES Q3! that the 
probability distribution of the local self energy, i?(E) [or equivalently of the diagonal element of the resolvent, Eq. (1) 
of the main text] is obtained in the large-TV limit from the equation: 


N 

= E* 

3=1 


ijGjj 


N 

E 

3=1 




E ^33 


( 20 ) 
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where = denotes the equality in distribution between random variables and E contains an infinitesimally small 
imaginary part i0 + . All correlations between the terms on the RHS of the previous relation can be neglected in the 
thermodynamic limit. This leads to a self-consistent equation on i?(E). For a full explanation of its derivation, see [1]. 
Here we just sketch how one can obtain an identity on its generating function. 

The first idea is that since the correlations between the terms on the RHS can be neglected, as it was shown in [2] , 
S.ji is a sum of a large number or heavy tailed i.i.d. variables and, hence, it’s a complex Levy random variable. We 
consider its generating function: 


dER(E)e- iXlS+iX2S * = e~ ih ij 


-ih\AX 1 G ij -X 2 G‘ j ) 


where X\ and X 2 are two real variables. The RHS can be computed in the following way: 

d h 


Y[ e - ih U x i G ^~ X2G W = (1 + I dER(E) 


2\h\ 1 +i* L 


— _ - V 2 ) 

3 U-E E-£* ) _ 


N 


Henceforth we neglect all the subleading (vanishing) terms in 1/N. This allows one to derive the identity: 


dER(E) e 


— iX\^2-\-iX2 S* _ 


= exp /r / dEi?(E) 


d h 


2 \h\ 


l+/i 


-nh 2 ( J£l _ X 2 \ 

UL V E-H E--E* ) — 


( 21 ) 


This is an implicit version of the self-consistent equation satisfied by R( E), which also defines R( S) as a complex 
Levy stable distribution [2]. 

We show in the following that this same result also follows directly from the saddle point equation (19) on p(<i>). 
By extremizing the action ( |18[ ) on /?($) and taking into account the normalization condition on p($), at leading order 
in N one finds: 

p(&) = exp (^E&£<f> + JdVC(&£V)p(*)\ . 


By plugging the expression (19) into the previous equation, one can perform the integral over 


dER(E) exp ^$ (1)t $ (1) (L;- E) - l$( 2 )t$( 2 )(£ _ E*)J = 
exp (^-E +fi J dsi? ( E ) J_ 


Ah 


-ih 2 




2\h\ 1 +f Jj 

This expression has to be valid for any and $( 2 )t$( 2 ), hence it defines a self-consistent on i?(E) which 

actually coincides with Eq. (21) established previously. This result shows that our super-symmetric formalism is in 
agreement with previous exact results: i?(E) is the complex Levy stable distribution obtained rigorously in [2]. 


DYSON BROWNIAN MOTION AND THE REGIME 1 < p < 2 

As we explained in the main text, we can show that the level statistics for 1 < p < 2 is the one of GOE matrices 
under the hypothesis that all Wigner matrices with the same heavy tails are characterized by the same level statistics 
(a very reasonable assumption). 

Our strategy consists first in modifying the distribution P(hij) of the matrix elements into (1 — e)P(hij) + 
eW/ fl W ( A r 1 / 11 h l: j) where W(x ) is a Gaussian distribution with unit variance. This does not alter the fat tails of 
the matrix elements and allows focusing on random matrices T-L e = (1 — e)"H + eW where T~L is a LM and W a very 
small GOE matrix whose elements have exactly the same scaling with N than the ones of TL. The level statistics of 
Hr— and, by the previous assumption, of H —can be obtained using the Dyson Brownian Motion. Since this is a very 
well known technique, we refer to the literature for an introduction jl4l 1 15] and just reproduce the main steps and 
ideas needed for our argument. 

Let’s denote A i(t) the eigenvalues of the matrix H t = (1 — t)H + tW. For t = 0 these coincides with the eigenvalues 
of H and for t = e with the ones of H e - The idea behind the Dyson Brownian Motion technique is to transform the 
interpolation t = 0 —► e in a stochastic process on the eigenvalues. The main trick is the property that a Gaussian 
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variable can be considered as the sum of two independent Gaussian variables. This allows one to describe the change 
from the Ht ensemble to the Ht+dt ensemble as an addition of a GOE matrix (with a suitable variance) and a rescaling. 
Taking an infinitesimal dt allows one to use perturbation theory: the final result that goes under the name of Dyson 
Brownian Motion is a stochastic evolution equation on the eigenvalues. The eigenvalues of H are simply the initial 
conditions for this process, whereas the ones of H e are the values obtained after a “time” t = e. 

The DBM admits a stationary distribution for the Xi(t)s which is simply the GOE distribution. The crucial question 
is then whether the A i(t) have enough “time” to equilibrate to their equilibrium (GOE) probability measure. What 
was conjectured already by Dyson and proved in great generality in recent years m is that the relaxation timescale 
to obtain local equilibration in the bulk of the spectrum to GOE statistics scales as N 1 ^ /N (with the scaling we 
have considered in this work). 

In consequence, for any /i £ (1,2), for any finite value of t. in particular t = e, the statistics of the eigenvalues in 
the bulk of the spectrum converges to the GOE one in the large -N limit. Using the assumption that all matrices with 
the same heavy tails are characterized by the same level statistics, we then find GOE level statistics for all matrices 
H e , in particular H = H e - q. 


NUMERICAL RESULTS FOR p <= (1,2) 


In this section we provide several numerical results obtained from exact diagonalizations of Levy Matrices in the 
range £ (1, 2), for several system sizes N = 2 m , with m from 8 to 14. As explained in the main text, the data are 
averaged over (at least) 2 22_m realization of the disorder. The energy spectrum is resolved in 64 small intervals v, 
centered around the energies E v = (\ n ) v . 

In fig. [4] we plot q tpp as a function of E„ for = 1.5 and for different system sizes, averaged over samples and 
eigenstates within each energy window. Although at high energies our numerical data are still quite far from full 
convergence, it is clear that q tpp evolves towards the GOE universal value, q^oE = 2/7r in all energy windows as N 
is increased. 

In fig. 1 we show the probability distribution of the gap ratio, n(r), for p, = 1.5, for different system sizes, and 
for four different values of the energy. In the case of Poisson statistics the probability distribution of the gap ratio 
is given by II(r) = 2/(1 + r) 2 , while the counterpart of II(r) corresponding to GOE statistics has been computed 
exactly in m■ Level repulsion in the GOE spectra manifests itself in the vanishing of the probability distribution 
at r = 0. As expected, at small enough energy the entire probability distribution II(r) is described by GOE. Finite 
size effects are stronger at higher energies. In particular for E ~ 13.08 (bottom-right panel of [5]), II(r) is still quite 
far from convergence even at the largest system size. Nevertheless it is clear that the distribution of the gap ratio is 
slowly evolving towards the GOE distribution as N is increased. 

Numerical results on the Inverse Participation Ratios (IPR) are coherent with previous findings. The IPR of 
the eigen-function n is defined as: T 2jrl = YliLi l( ? 'l n )| 4 - I n fig- |o] we plot the energy dependence of the exponent 
/3 = (lnT2 jn )j//hi N describing the scaling of the typical value of the IPR with the system size. At small enough 
energies we find /3 — 1, corresponding to the standard scaling of the IPR for fully delocalized states. As mentioned 
above, finite size effects are stronger at higher energies. We indeed observe that /? decreases as the energy grows for a 
fixed system size N. Nevertheless, at fixed energy, /I increases as the system size is increased and seems to approach 
the standard value 1 in all energy windows. This is confirmed by the numerical results on the the support set, recently 
introduced in DU as a tough measure wave-functions ergodicity. For an eigenvector n with sites ordered according to 

|(z|n)| > |(i + l|n)|, it is defined as the sets of sites i < Se^ such that J2i=i l(*l n )| 2 < 1 — e < J2i =1 +1 l(*l n )| 2 - The 
scaling of (S^) for N —» 00 and e arbitrary small but finite allows to discriminate between the extended and the 
localized regimes, as Se is TV-independent for localized wave-functions while it diverges for TV — > 00 for delocalized 
states. The exponent /?' = ln(S’e"^) 1/ / In TV, describing the scaling of the support set at large TV is also shown in fig. [fi] 
Its behavior is very similar to the one of /3 described above. However the support set is apparently a sharper measure 
of wave-functions ergodicity compared to the IPR, as the values of (3' are much closer to 1 in all energy windows. 

Similar results are obtained for fj, = 1.1, confirming that for fj, £ (1,2) all eigenstates of LMs are extended and the 
level statistics is described by GOE in the whole spectrum. Nevertheless finite size effect become stronger as [i is 
lowered and can be extremely important at high energies, where one needs to consider relatively large TV to observe 
full converges towards GOE. 
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FIG. 4. In (q typ /<1 qoe) as a function of the energy E for different system sizes for fj, = 1.5. 



FIG. 5. Probability distribution of the gap ratio for p = 1.5, for different system sizes, and for four different values of 
the energy. The Poisson and GOE counterparts of n(r) are also shown. Top-left panel: E = 0.024; The entire probability 
distribution is described by the GOE. Top-right panel: E = 5.53; II(r) converges to the GOE distribution for N large enough. 
Bottom-left panel: E = 9.48; II(r) evolves towards the GOE distribution as N is increased, although we are not able to observe 
full convergence. Bottom-right panel: E = 13.08; II(r) is very far from convergence even for the largest system size considered, 
although it slowly evolves towards the GOE distribution. 


NUMERICAL RESULTS FOR p <= (0,1) 


According to the expression ([3]) of the probability distribution of the entries of our LMs, each row (or column) of H 
has O(N) elements of and 0(1) elements of 0(1), ensuring a well-defined thermodynamic limit. However, 

the largest element of the whole matrix (which contains N 2 terms) is of order N 1 / /i . As a consequence, the range of 
variability of the matrix elements goes from O^N -1 ^) to which is, for large enough system sizes and for 

/i < 1, extremely broad. This could affect the numerical precision of our results. In order to overcome this issue, we 
have introduced a cut-off on large matrix elements scaling as ATV 1 ^, where A is a constant much larger than 1. Since 
we are only interested in the properties of LMs for energies of 0(1), the presence of such cut-off does not have any 
influence on our numerical results (provided that A is large enough). 

Furthermore, in order to exploit the “sparse-like” character of LMs, we introduce a cut-off 7 (very small but finite) 
on small matrix elements which allows to transform T~L in a sparse Erdos-Renyi random matrix constituted by the 
backbone of large entries 1 13] . This allows to simplify and speed-up the numerical calculations, since numerical routines 
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FIG. 6. Exponents fi (continuous lines) and /3' (dashed lines) describing the scaling with N of the typical value of the Inverse 
Participation Ratio and of the support set as a function of the energy for p = 1.5. 


for exact diagonalization are faster for sparse matrices. The probability distributions of the entry thus becomes: 

P ( N ,A \hij) = p { N s ( h ij) + (i ~p { n) °('y< \ h a\ < ^v 1/m a) 2 ^ 1+m , 

where = 1 - 1/(^7^) and C^’ A) = p/[y- p - A r ” 1 A _M ]. 

We have performed exact diagonalizations of such random matrices for several system sizes N = 2 m , with m from 
8 to 15, and for A = 2 15 . Data are averaged over (at least) 2 22_m realization of the disorder. The energy spectrum is 
resolved in 64 small intervals u, centered around the energies E v = (\ n )„. In order to make sure that the cut-off on 
small entries is small enough to reproduce the 7 —> 0 limit, we have considered different values of 7 (7 = 10~ 3 , 10 —4 , 
and 5 • 10~ 5 ) and checked that the data become independent of it (within our numerical accuracy). 

In the following we complement the discussion of the main text with more details and results obtained from exact 
diagonalizations for p £ (0,1). We will focus more specifically on p = 0.5. Similar results are found for p = 0.8 and 
p = 0.3, although finite size effects becomes bigger as p is decreased and the crossover region gets broader. 

In fig. 0 we plot q typ as a function of E v for p — 0.5 and for different system sizes, averaged over samples and 
eigenstates within each energy window. For small (resp. large) energies we recover, as expected, the universal values 
qp P = 2/7T (resp. cfp P —>• 0) corresponding to GOE (resp. Poisson) statistics. As mentioned in the main text, the 
curves corresponding to different values of N seem to cross much before the localization transition, which can be 
computed analytically and should occur at E * ~ 3.85. However, our numerical data on q typ are extremely clean and 
allow to observe that the crossing point is actually slowly drifting towards higher values of the energy (and most 
probably converging to E* in the thermodynamic limit). 

In fig. § we show the probability distribution of the gap ratio, n(r), for p = 0.5, for different system sizes, and 
for four different values of the energy. As expected, for small enough energies (e.g., E ~ 0.016, top-left panel) the 
entire probability distribution is described by GOE statistics. Conversely, for high enough energies (e.g., E ~ 7.68, 
bottom-right panel), in the localized regime, the data nicely approach the Poisson distribution n(r) = 2/(1 + r ) 2 — 
except for very small values of r where convergence is exponentially slow due to finite size effects. For moderately 
high energies (e.g., E = 1.25, top-right panel), n(r) evolves towards the GOE distribution as N is increased, although 
we are not able to observe full convergence for the largest system size. Finally, for energies in the crossover region 
(e.g., E = 2.28, bottom-left panel), one seems to observe that n(r) is described by a stationary (i.e., N- independent) 
and non-universal—neither GOE nor Poisson—distribution, as observed in |T] . Nevertheless, if one analyzes carefully 
the numerical data, focusing, for instance, on the behavior of n(r) at small r, one realizes that n(r) evolves in a 
non-monotonic way: for system sizes smaller than the crossover size, N < N rn ~ 1200 (see the inset of fig. 3 of 
the main text), it evolves towards the Poisson distribution, while for large system sizes, N > N rn , it commences to 
approach the GOE distribution. However, it is evident that if one ignored the existence of the crossover scale, based 
on the bottom-left panel of fig. [8] one would certainly conclude that for intermediate energies a new and non-universal 
“mixed” level statistics is found. 

In fig. [9] we plot the energy dependence of the exponents /3 = (In T 2>n )„/In N and /3' = ln(Se n ^)„/lnIV’ describing 
the scaling with the system size of the typical value of the IPR and of the average support set respectively, for p = 0.5. 
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FIG. 7. In (q typ /qao E ) as a function of the energy E for different system sizes for fj, = 0.5. 



r 


r 


FIG. 8. Probability distribution of the gap ratio for fi = 0.5, for different system sizes, and for four different values of the energy. 
The Poisson and GOE counterparts of II(r) are also shown. Top-left panel: E = 0.016; The entire probability distribution 
is described by the GOE. Top-right panel: E = 1.25; II(r) evolves towards the GOE distribution as N is increased, although 
we are not able to observe full convergence. Bottom-left panel: E = 2.28; II(r) seems to be described by a IV-independent 
non-universal distribution. Bottom-right panel: E = 7.68; II(r) converges to the Poisson distribution for large N. 


The behavior of /? and f3' is coherent with previous results, at least for sufficiently small and sufficiently large energies. 
More precisely, one observes that, at fixed N, /? and /?' decrease as the energy is increased. Nevertheless, at fixed and 
small enough energy, they both grow with N and seem to approach the standard value 1 for N —> oo. Conversely, 
at fixed and large enough energy, in the localized regime, /3 and j3' decrease to zero as the system size is increased, 
implying that (Y 2 , n )v, (5'c "'*)v —> cst. As mentioned above, the support set provides a more precise measure of wave- 
function ergodicity compared to the IPR. In particular, the exponent [3 is much smaller than one already very far from 
the localization transition. In the crossover region one should expect that /3 and /3 7 show a non-monotonic behavior 
as a function of N on the crossover scale N m (E). However, our numerical data are too noisy to capture this behavior. 
In fact, numerics based solely on the IPRs are inconclusive and could be undoubtedly misinterpreted, especially for 
intermediate energies within the crossover regime, since they are affected by strong finite size effects. 
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FIG. 9. Exponents fi (continuous lines) and /3' (dashed lines) describing the scaling with N of the typical value of the Inverse 
Participation Ratio and of the support set as a function of the energy for y = 0.5. 


NUMERICAL SOLUTION OF THE SELF-CONSISTENT EQUATION FOR Q(G) 


In this section we provide more details on the numerical solutions of the self-consistent equation on the probability 
distribution of the diagonal elements of the resolvent matrix. 

In order to device an accurate and efficient algorithm to compute Q(G) it is convenient to make use, again, of the 
“sparse-like” character of LMs m- As explained above, each row (or column) of % has 0(N) elements of 0(N 1 / /i ) 
and 0(1) elements of 0(1). We thus introduce a cut-off 7 that separates large matrix elements (\hij\ > 7 ) from small 
ones (| h t j \ < 7 ). The backbone of large entries constitutes thus a sparse (Erdos-Renyi) RM, with average connectivity 
c 7 = 2 N j P(h)dh = 7 -Al (13, IS]. Since the probability distribution of small matrix elements has now a finite 


variance ex 2 = 2 i/„ h 2 P{h) d h = /xy 2 ^/[N {2 — /z)], their contribution to the sum in Eq. (1) of the main text 

can be handled using the (classical) central limit theorem, yielding the following self-consistent equation for Q 7 (G) 
(in the limit 7 —>• 0 ): 


00 « k / 

Q 7 (G) = E P"i (^0 / | d/ii-P(/ii)] S I G 

k—0 ** 2=1 \ 


- E 


n 7 


*<G>+ £>?<?< 


( 22 ) 


i=0 


where p 7 (fc) = e C 7 c 7 //c! is the Poisson distribution of the connectivity. This equation can be efficiently solved 
using a population dynamics algorithm m We have used a population of 2 26 elements, and computed Q 7 (G) for 
7 = 10~ 3 ,10 -4 , 5 • 10” 5 , and extrapolated the results for 7 —> 0. 

In fig. [TO] we show the marginal probability distribution of In $sGu for several values of the imaginary regulator 
77 and for E = 20, deep in the GOE ergodic phase. Since the system is delocalized and the spectrum is absolutely 
continuous, Q/(ln3G) must have a non-singular limit as 77 —> 0 + . We indeed observe a stationary 77 -independent 
distribution for 77 sufficiently small (77 < 10 -6 ). As a consequence, (T 2 ) —1 0 for 77 — > 0 + . 

Conversely, in the localized phase the marginal probability distribution of the imaginary part of Gu has a singular 
behavior as 77 —► 0 + , as illustrated in fig. 11 Almost all values of ‘SsGa are of order 77, except extremely rare events— 
whose fraction vanishes as 77—described by heavy power-law tails with an exponent 1 + m and m = 1 / 2 . More 
precisely, Q/( 3 G) has a scaling form f(x/ 7 j)rj for x ~ 77, with Jf{y) d y = 1 , and fat tails 


Qii^G) ~ 


c? 7 


1 —m 


(3G) 


1+ra 


( 23 ) 


with c being a constant of 0 ( 1 ), and a cut-off for 3 G^ ~ I/77. Such tails gives a contribution of 0 ( 1 ) to the density of 
states, whereas the bulk part only yields a vanishing contribution. The marginal probability distribution of the real 
part of Gu (not shown) converges to a stationary distribution with power-law tails with a E-independent exponent 
1 + to = 2 . This implies that in the localized phase (T2) —> cst for rj —> 0 + . 


In fig. 


12 


that the 77 


we show the behavior of Qj(\n^G) in the crossover region. Since the system is delocalized, we know 
> 0 + limit exists and is non-singular. However, convergence to a stationary distribution is observed only 
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FIG. 10. Marginal probability distribution of ln9G for different values of the imaginary regulator 77 and for E = 1.25, showing 
convergence to a stationary ^-independent distribution for small enough 77 . 



FIG. 11. Marginal probability distribution of lnQG for different values of the imaginary regulator 77 and for E = 5.5. In the 
localized phase the limit 77 —> 0 + is singular: Almost all values of QGa are of order 77 , except extremely rare events described 
by heavy power-law tails with an exponent 1 + m = 3/2 whose coefficient vanishes as ^/rj. 


for extremely small values of the imaginary regulator, 77 < 10 “ 13 in this case. For 77 small enough but still larger 
than 10 —13 , one observes that, similarly to the localized regime, the marginal distribution of ^sGa displays “singular” 
power-law tails described by Q/(3G) ~ rj 1 2 3 4 5 ~ m /(3G) 1+m with an exponent 1/2 < 777 < 1, and a cut-off for ^sGa ~ I /77 
(the exponent is instead 771 ~ 1 for the marginal distribution of ReGjj). This implies that for large enough 77 the tails 
of Q/(lnSG) give a 0(1) contribution to the density of states, whereas the bulk part gives a contribution of 0( 77 ), as 
if the system was non-ergodic. 

In order to extract the crossover scale N' m (E), we measure the typical value of the imaginary part of G„;, SG-/ P = 
g(in <sGa ), over the stationary distribution on the loca lized phase and, following the argument presented in the main 


text, we define N' m (E) = 1/(SG*,- P p(i?)). In fig. 13 we plot In N' m (E) as a function of InN m (E), showing a linear 


relation between these two quantities. This implies that our argument allows to capture correctly the origin of finite 
size effects. 
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is achieved for 77 < 10 ~ 13 . 



l°g 2 N m 


FIG. 13. In N! m (E) as a function of In Nm(E). 


[ 6 ] R. Abou-Chacra, P.W. Anderson, and D.J. Thouless, J. Phys. C: Solid St. Phys. 6 , 1734 (1973). 

[7] V. Bapst, J. Math. Phys. 55, 092101 (2014). 

[ 8 ] K. Efetov, Supersymmetry in Disorder and Chaos, Cambridge Univ. Press 1997. 

[9] A.D. Mirlin and Y.V. Fyodorov, Nucl. Phys. B 366, 507 (1991). 

[10] E. Tarquini, G. Biroli, and M. Tarzia, in preparation. 

[11] A.D. Mirlin and Y.V. Fyodorov, J. Phys. I France 2, 1571 (1992). 

[12] Z. Burda, J. Jurkiewicz, M.A. Nowak, G. Papp, and I. Zahed, Phys. Rev. E 75, 051126 (2007). 

[13] F.L. Metz, I. Neri, and D. Bolle, Phys. Rev. E 82, 031135 (2010); JSTAT (2010) P01010. 

[14] For reviews see M.L. Mehta, Random Matrices, 3rd Edition, (Elsevier-Academic Press, 2004); S.N. Majumdar, Random 
Matrices, the Ulam Problem, Directed Polymers & Growth Models, and Sequence Matching, Les Houches lecture notes for 
the summer school on Complex Systems (Les Houches, 2006); T. Guhr, A. Mueller-Groeling, H.A. Weidenmueller, Phys. 
Rep. 299 189 (1998). 

[15] L. Erdos, A. Knowles, H.-T. Yau, J. Yin, Electronic J. of Probability 18, 1 (2013). 

[16] Y.Y. Atas, E. Bogomolny, O. Giraud, and P. Vivo, J. Phys. A: Math. Gen. 46, 355204 (2103). 

[17] A. De Luca, V.E. Kravtsov, B.L. Altshuler, and A. Scardicchio, arXiv: 1401.0019 

[18] K. Janzen, A. Engel, and M. Mezard, Europhys. Lett. 89, 67002 (2010). 

[19] M. Mezard and G. Parisi, Eur. Phys. J. B 20 (2001), 217. 








